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ABSTRACT 

Wc present a new method to retrieve molecular abundances and temperature profiles from exoplanct 
atmosphere photometry and spectroscopy. We run millions of ID atmosphere models in order to cover 
the large range of allowed parameter space. In order to run such a large number of models, we have 
developed a parametric pressure-temperature {P-T) profile coupled with linc-by-linc radiative transfer, 
hydrostatic equilibrium, and energy balance, along with prescriptions for non-equilibrium molecular 
composition and energy redistribution. The major difference from traditional ID radiative transfer 
models is the parametric P-T profile, which essentially means adopting energy balance only at the top 
of the atmosphere and not in each layer. We see the parametric P-T model as a parallel approach to the 
traditional exoplanet atmosphere models that rely on several free parameters to encompass unknown 
absorbers and energy redistribution. The parametric P-T profile captures the basic physical features 
of temperature structures in planetary atmospheres (including temperature inversions), and fits a wide 
range of published P-T profiles, including those of solar system planets. We apply our temperature 
and abundance retrieval method to the atmospheres of two transiting exoplanets, HD 189733b and 
HD 209458b, which have the best Spitzer and HST data available. For HD 189733b, we find efficient 
day-night redistribution of energy in the atmosphere, and molecular abundance constraints confirming 
the presence of H2O, CO, CH4 and CO2. For HD 209458b, we confirm and constrain the day-side 
thermal inversion in an average ID temperature profile. We also report independent detections of 
H2O, CO, CH4 and CO2 on the dayside of HD 209458b, based on six-channel Spitzer photometry. Wc 
report constraints for HD 189733b due to individual data sets separately; a few key observations are 
variable in different data sets at similar wavelengths. Moreover, a noticeably strong CO2 absorption in 
one data set is significantly weaker in another. We must, therefore, acknowledge the strong possibility 
that the atmosphere is variable, both in its energy redistribution state and in the chemical abundances. 

Subject headings: methods: numerical — planetary systems — planets and satellites: general - 
planets and satellites: individual (HD 209458b, HD 189733b) — radiative transfer 



1. INTRODUCTION 

Major advances in the detection and characterization 
of extrasolar planet atmospheres have been made during 
the last decade. Over a dozen hot Jupiter atmospheres 
have been observed by the Spitzer Space Telescope and a 
handful by the Hubble Space Telescope (HST). These ob- 
servations have changed the field of exoplanets, enabling 
for the first time studies of atmospheres of distant worlds. 

Observational highlights include the identification of 
molecules and atoms, and signatures of thermal inver- 
sions. HST observed water vapor and methane in trans- 
mission in HD 189733b (Swain et al. 2008b). Water 
vapor was also inferred in transmission at the 3.6 /im, 
5.8 /im, and 8 /im in HD 189733b (Tinetti et al. 2007) 
(but c.f. Desert et al. 2009, and also Beaulieu et al. 
2009). Water was inferred on the dayside of HD 189733b 
from the planet-star flux contrast in the 3.6 /im - 8.0 /im 
Spitzer broadband photometry (Barman 2008; Charbon- 
neau et al. 2008). Additionally, Grillmair et al (2008) 
reported a high signal-to-noise spectrum of the dayside of 
HD 189733b in the 5 /im - 14 /jm range, using the Spitzer 
Infrared Spectrograph (IRS), and reported detection of 
water. More recently, HST detected carbon dioxide in 
thermal emission in HD 189733b (Swain et al. 2009a), 
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in 1.4 nm. - 2.6 /im. Previously, sodium was detected in 
HD 209458b (Charbonneau et al. 2002) and HD 189733b 
(Redfield et al. 2008). Several observations of the plan- 
etary dayside have also been reported for HD 209458b, 
first observed by Spitzer at 24 /im (Deming et al. 2005 
and Seager et al. 2005). A major observational discovery 
was the finding of strong emission features in the broad- 
band photometry of HD 209458b at secondary eclipse 
(Knutson et al. 2008), implying a thermal inversion in 
the atmosphere (Burrows et al. 2007 & 2008). In tandem 
with our present work, water, methane and carbon diox- 
ide were independently inferred and recently reported on 
the dayside of HD 209458b, with ffSTNICMOS and five- 
channel Spitzer photometry (Swain et al. 2009b). Ad- 
ditionally, there are hints that variability of hot Jupiter 
thermal emission may be common (Knutson et al. 2007 
vs. Charbonneau et al. 2008; Charbonneau et al. 2008 
vs. Grillmair et al. 2008; Deming et al. 2005 vs. Dem- 
ing, private communication 2009). The discoveries of at- 
mospheric constituents and temperature structures mark 
the remarkable successes of exoplanet atmosphere mod- 
els in interpreting the observations. 

Despite the successes of model interpretation of data, 
major limitations of traditional self-consistent atmo- 
sphere models are becoming more and more apparent. 
Model successes include the detection of thermal inver- 
sions on the day-side of hot Jupiters (Burrows et al. 2007; 
Knutson et al. 2008), and subsequent classification of 
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systems on the same basis (Burrows et al. 2008; Fortney 
et al. 2008; Seager et al. 2008). The nature of the ab- 
sorbers causing the inversions, however, is not yet known 
(Knutson et al. 2008 & 2009a), forcing modelers to use 
an adhoc opacity source to explain the data. The intially 
favored opacity sources, TiO and VO may be unlikely to 
cause inversions of the observed magnitude (Spiegel et 
al. 2009). 

A second model limitation arises because ID radiative- 
convective equilibrium models cannot include the com- 
plex physics involved in hydrodynamic flows. Existing 
models use a parameterization of energy transfer from 
the day side to night side (e.g., Burrows et al. 2008), 
using parameters for the locations of energy sources and 
sinks, and for the amount of energy transported. 

One further example of atmospheric model limitations 
involves the treatment of chemical compositions. Self- 
consistent models generally calculate molecular mixing 
ratios based on chemical equilibrium along with the as- 
sumption of solar abundances, or variants thereof (e.g., 
Seager et al. 2005). However, hot Jupiter atmospheres 
should host manifestly non-equilibrium chemistry (e.g., 
Liang et al. 2003, Cooper & Showman, 2006), which ren- 
der the assumption of chemical equilibrium only a fidu- 
cial starting point. 

With parameters used to cover complicated or un- 
known physics or chemistry, existing "self-consistent" 
radiative-convective equilibrium models are no longer 
fully self consistent. Furthermore, given the computa- 
tional time of such models, only a few models are typ- 
ically published to interpret the data, often only barely 
fitting the data — a quantitative measure of fit is typically 
absent in the literature. In an ideal world, one would con- 
struct fully self-consistent 3D atmosphere models (e.g., 
Showman et al. 2009) that include hydrodynamic flow, 
radiative transfer, cloud microphysics, photochemistry 
and non-equilibrium chemistry, and run such models over 
all of parameter space anticipated to be valid. Such mod- 
els will remain idealizations until computer power im- 
proves tremendously, given that atmospheric circulation 
models can take weeks to months to run even with simple 
radiative transfer schemes. 

We are motivated to develop a data-interpretation 
framework for exoplanet atmospheres that enables us to 
run millions of models in order to constrain the full range 
of pressure-temperature (P-T) profiles and abundance 
ranges allowed by a given data set. In this work, we 
report a new approach to ID modeling of exoplanet at- 
mospheres. Taking a cue from parameterized physics al- 
ready being adopted by existing models, we go much fur- 
ther by parameterizing the P-T profile as guided by basic 
physical considerations and by properties of P-T profiles 
of planetary atmospheres in the literature. Indeed, it ap- 
pears at present, the only way to be able to run enough 
models to constrain the P-T structure and abundances 
is to use a parameterized P-T profile. The essential dif- 
ference between our new method and currently accepted 
radiative-convective equilibrium models is in the treat- 
ment of energy balance. We ensure global energy bal- 
ance at the top of the atmosphere, whereas conventional 
atmosphere models use layer-by-layer energy balance by 
way of radiative and convective equilibrium. Our new 
method can be used as a stand-alone model, or it can be 
used to identify the parameter space in which to run a 



reasonable number of traditional model atmospheres. We 
call our method "temperature and abundance retrieval" 
of exoplanet atmospheres. 

Atmospheric temperature and abundance retrieval is 
not new (e.g. Goody & Yung, 1989). Studies of plan- 
etary atmospheres in the solar system use temperature 
retrieval methods, but in the present context there is 
one major difference. Exquisite data for solar system 
planet atmospheres means a fiducial pressure - temper- 
ature profile can be derived. The temperature retrieval 
process for atmospheres of solar system planets therefore 
involves perturbing the fiducial temperature profile. For 
cxoplancts, where data are inadequate, there is no start- 
ing point to derive a fiducial model. Nevertheless, Swain 
et al. (2008b, 2009a) used previously published model 
P-T profiles of HD 189733b, and varied molecular abun- 
dances to report model fits to HST NICMOS spectra of 
HD 189733b. In another study, Sing et al. (2008) used 
a parametric temperature profile with a thermal inver- 
sion, and parametrized abundances, required specifically 
to explain HST STIS optical observations of HD 209458b. 
Our model, however, is completely general in the choice 
of P-T profiles, ranges over tens of thousands of profiles, 
and satisfies the physical constraints of global energy bal- 
ance and hydrostatic equilibrium. Our model reports, 
with contour plots, the quantitatively allowed ranges of 
P-T profiles and molecular abundances. In addition, we 
also report constraints on the albedo and day-night en- 
ergy redistribution, and on the effective temperature. 

In this paper, the P-T profiles are introduced in § 2, 
and the radiative transfer model and parameter space 
search strategy are described in § 3. Results of data 
interpretation with the model are presented in § 4 for 
HD 189733b, and in § 5 for HD 209458b. § 6 discusses the 
overall approach, and prospects for future work. And, in 
§ 7 we present a summary of our results and conclusions. 

2. TEMPERATURE - PRESSURE STRUCTURE 

Our goal is to introduce a new method of parametriz- 
ing one-dimensional (TD) models of exoplanetary at- 
mospheres. Our approach is to parametrize the 
pressure-temperature (P-T) profiles, instead of sepa- 
rately parametrizing each of the different physical pro- 
cesses that contribute to the P-T profiles. In what fol- 
lows, we will first review the basic physics behind the 
temperature structure of planetary atmospheres which 
serves as motivation for our P-T form. Next, we will de- 
scribe our parametric P-T profile. Finally, we compare 
our model P-T profiles with published profiles for some 
known planetary systems. 

2.1. Temperature Structure in Planetary Atmospheres 

The pressure-temperature profiles of planetary atmo- 
spheres are governed by basic physics. Here, we qualita- 
tively describe the physics and shape of a representative 
P-T profile starting from the bottom of the atmosphere. 
We focus on the generality that the temperature struc- 
ture at a given altitude depends on the opacity at that 
altitude, along with density and gravity. 

In the deepest layers of the planet atmosphere, con- 
vection is the dominant energy transport mechanism. 
The high pressure (cquivalcntly, high density) implies a 
high opacity, making energy transport by convection a 
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more efficient energy transport mechanism than radia- 
tion. For hot Jupiters, the daysidc radiative - convec- 
tive boundaries usually occur at the very bottom lay- 
ers of the atmospheres, at pressures higher than ~ 100 
bar. In the layers immediately above this boundary, the 
optical depth is still high enough that the diffusion ap- 
proximation of radiative transport holds. The diffusion 
approximation, and the incident stellar flux predomi- 
nantly absorbed higher up in the atmosphere, lead to 
a nearly isothermal temperature structure in the lay- 
ers immediately above the radiative-convective bound- 
ary in hot Jupiters. It is assumed here that the energy 
sources in the planet interior are weak compared to the 
stellar irradiation. In the solar system giant planet at- 
mospheres (and cooler exoplanet atmospheres as yet to 
be observed), the radiative - convective boundary occurs 
at lower pressures in the planet atmosphere than for hot 
Jupiters. 

Above the isothermal diffusion layer, i.e. at lower pres- 
sures, the atmosphere becomes optically thin and the 
diffusion approximation breaks down. These optically 
thin layers are at altitudes where thermal inversions can 
be formed, depending on the level of irradiation from the 
parent star and the presence of strong absorbing gases or 
solid particles. Thermal inversions, or "stratospheres" , 
are common to all solar system planets and have re- 
cently been determined to exist in several hot Jupiter 
atmospheres (Burrows et al. 2008). For the solar sys- 
tem giant planets, photochemical haze due to CH4 is the 
primary cause of stratospheres, and for Earth it is O3. 
For hot Jupiters, TiO and VO may be responsible for 
the thermal inversion (but see Spiegel et al. 2009), but 
the identification of the absorbers is still under debate. 
At still lower pressures, below P ~ 10~ 5 bar, the optical 
depths eventually become so low that the layers of the at- 
mosphere are transparent to the incoming and outgoing 
radiation. 

2.2. A Parametric P-T Profile 

We propose a parametric P-T profile, motivated by 
physical principles, solar system planet P-T profiles, and 
ID "self-consistent" exoplanet P-T profiles generated 
from model atmosphere calculations reported in the liter- 
ature. Our "synthetic" P-T profile encapsulates strato- 
spheres, along with the low and high pressure regimes 
of the atmosphere. We will refer to our synthetic P-T 
profile as the "parametric P-T profile" . 

The atmospheric altitudes of interest are those where 
the spectrum is formed. Nominally, we consider this 
range to be between 10 -5 bar < P < 100 bar, and we 
will refer to this pressure range as the "atmosphere" in 
the rest of the paper. At pressures lower than 10 -5 bar, 
the atmosphere layers are nearly transparent to the in- 
coming and outgoing radiation at visible and infrared 
wavelengths. And, above pressures of ~ 100 bar (or even 
above pressures of ~ 10 bar), the optical depth is too 
high for any radiation to escape without being repro- 
cessed. The corresponding layers do not contribute sig- 
nificantly to the emergent spectrum. Figure [T] shows a 
schematic parametric P-T profile. In our model, we di- 
vide the atmosphere into three representative zones or 
"Layers" as shown in Figure [1] The upper-most layer, 
Layer 1, in our model profile is a "mesosphere" with no 
thermal inversions. The middle layer, Layer 2, represents 
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Fig. 1. — The parametric pressure-temperature (P-T) profile. 
In the general form, the profile includes a thermal inversion layer 
(Layer 2) and has six free parameters, see . For a profile with 
thermal inversion, P2 > Pi, and for one without a thermal in- 
version, Pi > P2. An isothermal profile is assumed for pressures 
above P3 (Layer 3). Alternatively, for cooler atmospheres with no 
isothermal layer, Layer 2 co uld extend to deeper layers and Layer 
3 could be absent (see § 12.311 . 

the region where a thermal inversion (a "stratosphere") 
is possible. And, the bottom-most layer, Layer 3, is the 
regime where a high optical depth leads to an isothermal 
temperature structure. Layer 3 is used with hot Jupiters 
in mind; for cooler atmospheres this layer can be absent, 
with Layer 2 extending to deeper layers. 

Our proposed model for the P-T structure in Layers 1 
and 2 is a generalized exponential profile of the form: 

where, P is the pressure in bars, T is the temperature in 
K, and P , T Q , a and (3 are free parameters. For Layer 3, 
the model profile is given by T = T3, where T3 is a free 
parameter. 

Thus, our parametric P-T profile is given by: 

P <P<Pi : P = P e ai< - T - T °^ 1 Layerl 
P 1 <P<P 3 : p = p 2e »2(T-T 2 )^ Layer2 (2) 
P > P 3 : T = T 3 Layer3 

In this work, we empirically set /3i = fa = 0.5 (see 
i j2.3[) . The space of P-T profiles is spanned in all gener- 
ality with the remaining parameters, rendering /3 redun- 
dant as a free parameter. Then, the model profile in ([2|) 
has nine parameters, namely, Pq, To, a%, Pi, Pi, T2, a%, 
P3, and T3. Two of the parameters can be eliminated 
based on the two constraints of continuity at the two 
layer boundaries, i.e., Layers 1-2 and Layers 2-3. And, 
in the present work, we set Pq = 10 -5 bar, i.e., at the top 
of our model atmosphere. Thus, our parametric profile 
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Fig. 2. — Comparison of the parametric P-T profile with previously pub lished profiles. In each case, the dotted line is the published 
P-T profile and the solid line is a fit with the parametric profile (see § 12.311 . The P-T profiles of solar-system planets (upper- left panel) 
were obtained from the NASA Planetary Data System (see text for references). In the lower-right panel, "B08", "F06" and "K09" refer to 
Burrows et al. (2008), Fortney et al. (2006) and Knutson et al. (2009a), respectively. 



in its complete generality has six free parameters. 

Our P-T profile consists of fOO layers in the pressure 
range of 10 -5 — 100 bar, encompassing the three zonal 
"Layers" described above. The 100 layers are uniformly 
spaced in log(P). For a given pressure, the temperature 
in that layer is determined from equation using the 
form T = T(P). The kinks at the layer boundaries are 
removed by averaging the profile with a box-car of 10 
layers in width. 

2.3. Comparison with known P-T Profiles 



The parametric P-T profile described in § 12.21 fits the 
actual temperature structures of a wide variety of plane- 
tary atmospheres. Figure [2] shows the comparison of our 
model P-T profile with published P-T profiles of several 
solar system planets and hot Jupiters. For each case, 
the published profile was fitted with our six-parameter 
P-T profile using a Levenberg-Marquardt fitting proce- 
dure (Levenberg, 1944; Marquardt, 1963). We were able 
to find best fits to all the published profiles with the /3 
parameter fixed at 0.5; indicating that f3 is a redundant 
parameter. Therefore, we set fti = (3% = 0.5 for all the 
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P-T profiles in this work. In addition, it must be noted 
that, for the parametric P-T profile to have a thermal 
inversion, the condition P2 > Pi must be satisfied, and 
for one without a thermal inversion, Pi > P2. 

For the solar system planets (top left panel of Fig- 
ure [5]), published profiles show detailed temperature 
structures obtained via several direct measurements and 
high-resolution temperature retrieval methods (Lindal 
et al. 1981; Lindal et al. 1985; Lindal et al. 1987; 
Lindal et al. 1990; U.S. Standard Atm osphere 1976; 
|http:/ /atmos. nmsu.edu/planetary _datasets/ [). For hot 
Jupiters, on the other hand, observations are limited. 
Consequently, the published profiles for hot Jupiters 
(Figure 0) are obtained from self-consistent 1-D mod- 
els reported in the literature attempting to explain the 
available data on each system. As is evident from Fig- 
ure [H the parametric P-T profile provides good fits to 
all the known profiles. We emphasize that, while the fits 
alone are not a test of our retrieval method, they show 
that our parametric P-T profile is general enough to fit 
the wide range of planetary atmospheres represented in 
Figure El 

We also tested the emergent spectrum obtained with 
our parametric P-T profile by comparing it with a self- 
consistent model. In this regard, we used a P-T profile 
from a self-consistent atmosphere model (Scagcr et al. 
2005) to generate an emergent spectrum for HD 209458b. 
The emergent spectrum calculated by our model (for the 
same molecular composition) agrees with that of the self- 
consistent model. 

The fundamental significance of our parametric P-T 
profile is the ability to fit a disparate set of planetary 
atmosphere structures, with very different atmospheric 
conditions (Figure [2]). In particular, the published hot 
Jupiter P-T profiles were obtained from several different 
modeling schemes. The different planetary atmospheres 
all conform to a basic mathematical model because of 
the common physics underlying the P-T profiles. We 
anticipate that our parametric P-T profile will open up 
a new avenue in retrieving P-T structure of exoplanctary 
atmospheres by enabling an efficient exploration of the 
temperature structures and compositions allowed by the 
data. 

An important point concerns the adaptability of the 
model proposed in this work. The introduction of the 
isothermal Layer 3 in ([2]) is motivated by the fact that 
hot Jupiters, which are the prime focus of this work, 
have convective-radiative boundaries that are deep in the 
atmosphere, and the presence of an isothermal Layer 3 
is physically plausible (as is evident from self-consistent 
calculations in literature). However, the model can 
be adapted to cases of cooler atmospheres where the 
convective-radiative boundary occurs at lower pressures 
than for hot Jupiters, and where there may not be an 
isothermal layer. The appropriate model P-T profile in 
such cases would be one with only the two upper layers 
(Layer 1 and Layer 2). As can be seen from Figure [2 
model fits to the solar system planets belong to the cat- 
egory of models with only Layers 1 and 2. 

3. MODEL ATMOSPHERE 

The goal of our model atmosphere is to be able to gen- 
erate atmospheric spectra for a given planet over wide re- 
gions of parameter space and determine contours of fits 



with data. To meet this goal, our model uses efficient 
paramctrization of the pressure-temperature structure 
and chemical composition of a planetary atmosphere. 
Our model atmosphere includes a ID parametric P-T 
profile and parametric molecular abundances, coupled 
with line-by-line radiative transfer, hydrostatic equilib- 
rium, and the requirement of energy balance at the top 
of the atmosphere. Given a set of parameter values, the 
output is the emergent spectrum, and a goodness-of-fit 
to a given data set. 

The major difference of our model from traditional at- 
mosphere models is in the treatment of energy balance. 
Our model requires energy balance at the top of the at- 
mosphere, instead of an iterative scheme to ensure layer- 
by-laycr radiative (or radiative + convective) equilibrium 
as is done in conventional models. For a given set of 
model parameter values, we require that the net energy 
output at the top of the atmosphere (obtained by inte- 
grating the emergent spectrum) balances the net energy 
input due to the incident stellar flux (see ^3.3[) . Mod- 
els which do not satisfy this requirement are discarded. 
By running a large number of (~ 10 7 ) models in the pa- 
rameter space, and discarding those that do not satisfy 
the above requirement, we are left with a population of 
models that satisfy energy balance. 

Our temperature and abundance retrieval method also 
differs from other retrieval methods (Swain et al. 2008b, 
2009a; Sing et al. 2008). Our work differs fundamen- 
tally by: using a parametric P-T profile, exploring a 
wide range in parameter space (~ 10 4 P-T profiles inde- 
pendent of existing models, and ~ 10 models including 
abundances), maintaining hydrostatic equilibrium and 
global energy balance, and placing constraints on the 
day-night energy redistribution. We, therefore, consider 
our method to be a new step towards atmospheric tem- 
perature and abundance retrieval for extrasolar planets. 

3.1. Radiative Transfer and Hydrostatic Equilibrium 

We calculate emergent spectra using a one-dimensional 
line-by-line radiative transfer code, assuming LTE. Our 
model atmosphere consists of 100 layers in the pressure 
range between 10 -5 bar and 100 bar, and uniformly 
spaced in log(P). The key aspect of our model is the 
parametrization of the pressure-temperature profile and 
the chemical composition. In each layer, the parametric 
P-T profile (described in § 12. 2|) determines the temper- 
ature and pressure. And, the density is given by the 
ideal gas equation of state. The radial distance is de- 
termined by the requirement of hydrostatic equilibrium. 
And, the equation of radiative transfer, assuming LTE 
and no scattering, is then solved along with the adopted 
molecular abundances (sec § l3.2l bclow) to determine the 
emergent flux. The geometry and the radiative trans- 
fer formulation are modified appropriately for calculating 
transmission spectra. 

Our model directly computes the thermal emission 
spectrum and indirectly accounts for stellar heating and 
scattered radiation. In our model, we assume that the 
effects of heating by stellar radiation are included in the 
allowed "shapes" of the parametric P-T profile. As de- 
scribed in § 12.31 our parametric P-T profile fits the P-T 
profiles obtained from self-consistent models that directly 
include absorbed stellar radiation, including profiles with 
thermal inversions. 
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One reason the approach of indirectly including stel- 
lar irradiation is valid is that the visible and infrared 
opacities arc largely decoupled (Marlcy ct al. 2002). For 
example. Na, K, TiO, VO can have strong absorption 
at visible wavelengths in brown dwarfs and exoplanets, 
yet are mostly negligible absorbers in the infrared. H2O 
dominates infrared absorption, and CH4 and CO are also 
spectroscopically active in the IR. Another reason the 
indirect account for stellar irradiation via a parametric 
P-T profile is justified is the uncertain chemistry (via un- 
known optical absorbers that in some cases cause ther- 
mal inversions) and physics that is hard to account for 
directly in a self-consistent ID model. We note that it is 
not known whether the common assumption of radiative 
equilibrium in self-consistent models, while one possible 
way of including stellar irradiation, holds for hot Jupiter 
atmospheres. 

We include molecules of hydrogen (H 2 ), water vapor 
(H2O), carbon monoxide (CO), carbon dioxide (CO2), 
methane (CH4), and ammonia (NH3). Our H2O, CH4, 
CO and NH3 molecular line data are from Freedman et 
al. 2008, and references therein. Our CO2 data are from 
Freedman (personal communication, 2009) and Roth- 
man et al. (2005). And, we obtain the H 2 -H 2 collision- 
induced opacities from Borysow ct al. (1997), and Bo- 
rysow (2002). 

3.2. Molecular Abundances 

Molecular abundances in cxoplanct atmospheres may 
depart from chemical equilibrium (Liang et al. 2003; 
Cooper & Showman, 2006). Our code, therefore, has 
parametric prescriptions to allow for non-equilibrium 
quantities of each molecule considered. We begin by cal- 
culating a fiducial concentration of each molecule in a 
given layer, based on the assumption of chemical equi- 
librium. For chemical equilibrium we use the analytic ex- 
pressions in Burrows and Sharp (1999), and assume solar 
abundances for the elements. The molecules we include 
arc H 2 0, CO, CH 4 , C0 2 and NH 3 . The molecule C0 2 
is thought to originate from photochemistry and no sim- 
ple expression for equilibrium composition is available; 
we therefore choose an arbitrary fiducial concentration. 
We denote the concentration of each molecule by its mix- 
ing ratio, i.e number fraction with respect to molecular 
hydrogen 

The concentration of each molecule in a given layer is 
calculated by multiplying the corresponding fiducial con- 
centration by a constant factor specific to each molecule. 
The constant factor corresponding B to each molecule is 
a parameter of the model. Therefore, corresponding to 
the four prominent molecules, H2O, CO, CH4 and CO2, 
in our model, we have four parameters: /h 2 Oj /cOj /ch 4 
and /co 2 - For instance, /h 2 o is the ratio of the con- 
centration of H 2 to the fiducial concentration of H 2 
in each layer. Thus, although the concentration of H2O 
is different in each layer, it maintains a constant ratio 
(/h 2 o) with the fiducial concentration in that layer. By 
varying over a range in these four parameters, we are es- 
sentially varying over a wide range of molecular mixing 
ratios and elemental abundances. 

For NH3, we restrict the concentration to the fiducial 
equilibrium concentration, i.e /nh 3 = 1- This is because, 
NH3 has limited and weak spectral features in the Spitzer 
and HST bands under consideration in this work. Addi- 



tionally, NH3 is not expected to be abundant in enhanced 
quantities at such high temperatures as are seen in hot 
Jupiters. 

Clouds are not included in our model because HD 
209458b and HD 189733b likely do not have clouds that 
scatter or emit at thermal infrared wavelengths. Several 
reasons have been proposed in literature justifying the 
use of cloud-free models for these planets: finding weak 
effects of clouds on the P-T profile (Fortney et al. 2006 & 
2008); fast sedimentation rates in radiative atmospheres 
(Barman et al. 2005); clouds forming too deep in the 
planetary atmosphere to affect the emergent spectrum 
(Fortney et al. 2006); spectral features in transmission 
dominated by molecular absorption over cloud particles 
(Showman et al. 2009 and references therein); haze ab- 
sorption at blue wavelengths on HD 189733b drops off 
rapidly with increasing wavelength and should be negli- 
gible at IR wavelengths (Pont et al. 2008). That said, 
clouds can be included in our model as another free pa- 
rameter, via a wavelength-dependent opacity, and we 
plan to make this extension in the future. 

3.3. Energy Balance 

A planetary atmosphere model must satisfy the fun- 
damental constraint of energy conservation (here called 
energy balance). We consider a ID plane-parallel at- 
mosphere with normal incidence. While using a plane- 
parallel model, it is imperative to weight the incident 
stellar flux appropriately in order to represent an aver- 
age day side atmosphere. 

The average stellar flux incident on the planetary day 
side is given by F s = fF* where, F± is the wavelength in- 
tegrated stellar flux at the sub-stellar point of the planet, 
/is a geometric factor by which the stellar flux at the sub- 
stcllar point must be weighted so as to represent an av- 
erage ID plane-parallel incidence. A value of / = | rep- 
resents uniform distribution of stellar flux on the planet 
day side; it comes from x (n R^) / '(2tt R^) , i.e., the ratio 
of the planetary surface projected perpendicular to the 
incident stellar flux to the actual area of the planetary 
surface receiving the flux. However, it has been shown 
that just before secondary eclipse the daysidc average 
flux should be biased towards a higher value of / = | 
(Burrows et al. 2008). Therefore, in the present work, 
we adopt a value of / = | . 

In the context of our ID atmosphere, the constraint 
of energy balance means that the wavelength integrated 
emergent flux at the top of the planetary atmosphere 
matches the wavelength integrated incident stellar flux, 
after accounting for the Bond albedo (As) and possible 
redistribution of energy onto the night side (f r )- Our 
prescription for the day-night redistribution is similar to 
that of the P n prescription of Burrows et al. (2006). We 
define f r as the fraction of input stellar flux that is redis- 
tributed to the night side. The input stellar flux is given 
by Fi n = (1 — Ab)F s , i.e., the incident stellar flux (F s ), 
less the reflected flux (AbF s ). Therefore, the energy ad- 
vected to the night side is given by f r Fi n . And hence, 
the amount of flux absorbed on the dayside is given by 
Fn.day = (l-/r)-Fin, which equals (1 - A B ) (1 - f r )F s , It 
is assumed that the intrinsic flux from the planet interior 
is negligible compared to the incident stellar irradiation. 

Thus, the energy balance requirement on the day side 
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spectrum is given as: 

F p = (l-A B )(l-f r )F Sl (3) 

where, F p is the wavelength integrated emergent flux cal- 
culated from the model spectrum. We use a Kurucz 
model for the stellar spectrum (Castelli & Kurucz. 2004). 
We adopt the stellar and planetary parameters from Tor- 
res et al (2008). 

We use energy balance, equation ([3]), in two ways. The 
first, as described above, is to ensure that our model 
atmosphere satisfies energy conservation. The constraint 
on the model comes from the fact that Ab and f r are 
bounded in the range [0,1]. Thus, models for which r/ = 
(1 — Ag)(l — f r ) is greater than unity are discarded on 
grounds of energy balance. Secondly, for models which 
satisfy energy balance, equation ([3]) gives us an estimate 
of r\ for a given model. If an Ab is assumed, this gives a 
constraint on f r , and vice- versa. 

3.4. Parameter Space 

The parameter space of our model consists of N = 
Nprofiie + N mo i ec free parameters, where N pro f Ue is the 
number of parameters in the P-T profile, and N mo i ec 
is the number of parameters pertaining to the chemi- 
cal composition. As described in § 2, our parametric 
P-T profile has N pro fu e = 6, the six parameters being 
To, ax, P±, Pi, ot2, and P3. And, corresponding to the 
four prominent molecules, H 2 0, CO, CH4 and CO2, in 
our model, we have N mo i ec = 4 parameters: /h 2 Oi /cOi 
/ch 4 and /co 2 ( see 33. 2|) . Thus, our model has ten free 
parameters. 

Given a set of observations of a planet, for example 
a day-side spectrum, our goal is to determine the re- 
gions of parameter space constrained by the data, and to 
identify the degeneracies therein. We investigate models 
with and without thermal inversions. To this effect, we 
run a large number of models on a grid in the param- 
eter space of each scenario. Each parameter in either 
scenario has a finite range of values to explore, and we 
have empirically determined the appropriate step size in 
each parameter. A grid of reasonable resolution in this 
ten-parameter space can approach a grid-size that is com- 
putationally prohibitive. 

Even before running the grid of models, however, the 
parameter space of allowed models can be constrained 
to some extent based on the data and energy balance. 
Firstly, the model planet spectrum is bounded by two 
blackbody spectra, corresponding to the lowest and high- 
est temperatures in the atmosphere. Correspondingly, 
the maximum and minimum temperatures, T max and 
T m in, of the P-T profile are partly constrained by the 
maximum and minimum blackbody temperatures admis- 
sible by the data, and allowing a margin accounting for 
the fact that observations are available only at limited 
wavelengths. Secondly, an upper limit can be placed on 
the temperature at the base of our model atmosphere (T3 
in the P-T profile) based on energy balance. The max- 
imum of T3 is given by the effective temperature of the 
planet day-side assuming zero albedo and zero redistri- 
bution of the energy to the night-side, and allowing for 
a factor close to unity (~ 1 — 1.5) to take into account 
the changes in the blackbody flux due to line features. 
Therefore, given the data, the space of P-T profiles can 



be constrained to some extent even before calculating the 
spectra, using the above conditions. 

The final grid is decided by running several coarse 
grids, and obtaining an empirical understanding of the 
parameter space. At the end, a nominal grid in ten pa- 
rameters, with or without inversions in the P-T profiles, 
consists of ~ 10 7 models. We thus computed a total of 
~ 10 7 models for each planet under consideration. We 
ran the models on ~ 100 2-GHz processors on a Beowulf 
cluster, at the MIT Kavli Institute for Astrophysics and 
Space Research, Cambridge, MA. 

3.5. Quantitative Measure of Error 

Broad-band observations of exoplanet atmospheres are 
currently limited to the six channels of broadband in- 
frared photometry using Spitzer. Additionally, a few 
narrow band spectrophotometry have also been reported 
using the HST NICMOS G206 grism (1.4 - 2.6 fim), and 
the Spitzer IRS spectrograph (5-14 ^m). The number of 
broadband observations available in a single data set are 
typically smaller than the number of parameters in our 
model. Also, given that the parameters in our model are 
highly correlated and the degeneracies not well under- 
stood, the nature of existing data does not allow spectral 
fitting in the formal sense. 

Consequently, our approach in this work is to com- 
pute a large number of models on a pre-defined grid in 
the parameter space. For each model, we calculate a 
"goodness-of-fit" with the data using a statistic defined 
by a weighted mean squared error, given by: 

^2 1 \ 1 J fi,model fi,obs \ 

Here, fi is the flux observable and N b s is the number of 
observed data points. fi. b s and fi, m odei are the observed 
fi and the corresponding model fi , respectively. For sec- 
ondary eclipse spectra, fi is the planet-star flux contrast, 
whereas, for transmission spectra, fi is the transmis- 
sion. 0~i^ o bs is the 1 a measurement uncertainty in the 
observation. In order to obtain the model fluxes in the 
same wavelength bins as the observation, we integrate 
the model spectra with the transmission functions of the 
instruments. Where spectrophotometry is available, we 
additionally convolve our model spectra with the instru- 
ment point-spread function. 

The goodness-of-fit (£ 2 ), as defined in (j4|, is similar in 
formulation to the conventional definition of the reduced 
X 2 statistic, if N obs is replaced by the number of degrees 
of freedom. In the current context, we refrain from us- 
ing the x 2 statistic to asses our model fits because our 
number of model parameters (N) are typically more than 
the number of available broadband data points, leading 
to negative degrees of freedom. A value of x 2 with nega- 
tive degrees of freedom lacks meaning in that we cannot 
relate it to a confidence level. Therefore, in using £ 2 as 
our statistic of choice, we are essentially calculating x 2 
per data point as opposed to x 2 per degree of freedom. 
We discuss this further in § 16.21 

4. RESULTS: HD 189733B 
4.1. Secondary Eclipse 
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HD 189733b has the best data set for any hot Jupiter 
known to date, including both signal-to-noise and wave- 
length coverage. Several remarkable observations of 
the planet day side have been made by determining 
the planet-star flux ratios just before secondary eclipse. 
Grillmair et al. (2008) reported an infrared spectrum of 
the dayside in the 5 /xm - 14 /im range, obtained with the 
Spitzer IRS instrument. Charbonneau et al. (2008) pre- 
sented broadband photometry in the six Spitzer channels 
(3.6, 4.5, 5.8, 8, 16, and 24 /xmj3. Separate photomet- 
ric observations in the Spitzer 8 /xm and 24 (im channels 
were also reported by Knutson et al. (2007) and (2009b), 
respectively. Additionally, an upper limit was placed on 
the 2.2 /xm flux constrast by Barnes et al (2007). More 
recently, Swain et al. (2009a) observed the planet day- 
side with HST NICMOS spectrophotometry in the 1.5 
/xm - 2.5 /xm range. 

Together, these data sets place important constraints 
on the dayside atmosphere of the planet. The high S/N 
IRS spectrum reported by Grillmair et al (2008), ob- 
tained with 120 hours of integration time, is the best 
exoplanet spectrum known. This spectrum forms the 
primary data set for our present study of HD 189733b. 
In addition, we also use the photometric observations of 
Charbonneau et al. (2008), and the spectrophotometric 
observations of Swain et al. (2009a) to place complemen- 
tary constraints 0. We use the two independent obser- 
vations of Knutson et al. (2007 & 2009b), at 8 /xm and 
24 /xm, only as a guide for the models, and exclude them 
from our fits. This is because we place constraints from 
each individual data set separately, for reasons explained 
below, and we cannot place any constraints from a single 
data point from each of Knutson et al. (2007 & 2009b). 
However, it must be noted that the 8 //m and 24 /xm 
values of Knutson et al. (2007 & 2009b) agree with the 
IRS spectrum and the Charbonneau et al. (2008) obser- 
vations, respectively, which are being used in the current 
work. We present the data in Figure along with model 
fits. 

We analyze each data set separately. The reason we use 
a separate treatment is that some individual observations 
at similar wavelengths differ from each other. Specifi- 
cally, the 8 /xm Spitzer IRAC data point from Charbon- 
neau et al (2008) is ~ 2a above the IRS spectrum at 
the same wavelength (Grillmair et al. 2008). Also, 8 
/im IRAC observations taken at separate times differ by 
over 2 a (Charbonneau et al. 2008). From a model in- 
terpretation viewpoint, a 2a spread means that almost 
no constraints on molecular abundances can be drawn 
(based on the reported uncertainities in the data). In 
addition, the very strong evidence of CO2 that is visibly 
obvious around 2 /xm in the HST/NICMOS data set is 
not so strong in the 4.5 /xm and 16 /xm Spitzer channels. 
In this work, we are unable to find models that fit all 
the three data sets of HD 189733b simultaneously within 
the £ 2 = 2 level. The data and model results below may 
well indicate atmospheric variability in this planet. See 
§ 16.31 for further discussion about possible atmospheric 

3 The 16 /xm value of Charbonneau et al. (2008) was obtained 
by reanalyzing the observations of Deming et al. (2006), and has 
been found to be consistent with the latter. 

4 We do not use those binned data points of Swain et al. (2009a) 
that are either at the edges of the chip or are non-detections at the 
3-<T level. These include points with A < 1.65 /xm and A > 2.4 /jm 



variability. 

In what follows, we describe the constraints on the 
atmospheric properties of HD 189733b placed by each 
data set. Ideally, we would like to report constraints at 
the £ 2 = 1 surface, i.e. only those models which, on 
an average, fit the observations to within the la error 
bars. However, we adopt a more conservative approach 
and lay emphasis on the results at the £ 2 = 2 surface 
as well, allowing for models fitting the observations to 
within ~ 1.4cr on an average. The latter choice is some- 
what arbitrary, and the informed reader could choose to 
consider other £ 2 surfaces; we show contours of £ 2 in the 
different atmospheric properties. 

4.1.1. Molecular Abundances 

The constraints on the abundances due to the differ- 
ent data sets are shown in Figure |U The figure shows 
contours of £ 2 in the parameter space of our model at- 
mosphere. The contours show minimum £ 2 surfaces. At 
each point in a given two-parameter space, several val- 
ues of £ 2 are possible because of degeneracies with other 
parameters in the parameter space. In other words, the 
same point in a pair of parameters can have multiple 
values of £ 2 corresponding to evaluations at various val- 
ues in the remaining parameters. For the contours in 
Figure [U we consider the minimum possible value of £ 2 
at each point in a two-parameter space. Consequently, 
each colored surface shows the region of parameter space 
which allows a minimum possible £ 2 corresponding to 
that color. The mixing ratios of species are shown as 
ratios by number. 

Spitzer IRS spectrum: We begin with the IRS spec- 
trum (Grillmair et al. 2008) because this data set has 
the broadest wavelength coverage and the best S/N. Also, 
the IRS data set has more data points that the number 
of parameters in our model. Our most significant new 
finding is the existence of models fitting the IRS data set 
that are totally consistent with efficient redistribution 
of energy from the planet dayside to the night side (see 
Figure S] for the range of allowed values). This addresses 
a significant problem in previously published models for 
dayside spectra (Grillmair et al. 2008) not being able 
to explain the efficient redistribution found in the phase 
curves of Knutson et al. (2007 & 2009b) (See Charbon- 
neau et al. 2008, and Grillmair et al. 2008; but also see 
Barman 2008). With the IRS data, we also confirm the 
absence of a thermal inversion in HD 189733b. 

Despite the high S/N, the IRS spectrum allows only 
weak constraints on the molecular abundances. The rea- 
son being that the IRS spectral range has molecular ab- 
sorption features of only CH 4 and H 2 0; even for these 
molecules, their features at other wavelengths remain 
completely unconstrained. The significant features are 
of CH4 (7.6 /xm), and H2O (6.3 /im). For example, the 
£ 2 = 2 surface constrains the mixing ratio of water vapor 
to 10" 6 - 0.1, and that of CH 4 to be less than 10~ 2 (but 
c.f §E3J. 

Spitzer broadband photometry: Our model fits to the 
Spitzer broadband data (Charbonneau et al. 2008) re- 
quire less efficient circulation of incident stellar energy 
compared to the best model fits to the Spitzer IRS spec- 
trum. Specifically, the value for redistribution in our 
framework is f r < 0.25, and this agrees with previous 
results (Barman 2008, Charbonneau et al. 2008) which 
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Fig. 3. — Secondary eclipse spectra for HD 189733b. The black squares (filled with red) with error bars show the data from Spitzer 
photometry in six channels, at 3.6, 4.5, 5.8, 8, 16, and 24 /im, reported by Charbonncau et al. 2008 (the 16 (im point is from reanalysis of 
Doming ct al. 2006, consistent with the latter). The black circles with error bars in the main panel show the IRS data from Grillmair et al. 
2008, and those in the inset show the HST NICMOS data from Swain et al. 2009a (see § I4.1l for details). The upper-limit at 2.2 /im shows 
the 1-(T constraint from Barnes et al. (2007). In the main panel, the red, orange, and green spectra show model spectra that fit the IRS 
spectrum with § 2 in the range of 1.0 - 2.0, 2.0 - 3.0, and 3.0 - 4.0, respectively. Ten spectra from each category are shown. For guidance, 
a single best fit model spectrum for this data is shown in blue. The gray spectra in the main panel show ten of the models that fit the 
Spitzer broadband photometry with £ 2 < 2.0. For guidance, a single best fit model spectrum for this data is shown in brown (with the 
bandpass averaged points shown in light-blue). The inset shows ten of the model spectra that fit the NICMOS data; with colors following 
those in the main panel. The blue and brown models in the inset are one each best-fit to the IRS data and photometric data, respectively, 
showing that the best-fit models for the IRS data and the broadband photometry do not match the NICMOS data. 



find that model fits to the broadband data require inef- 
ficient day-night redistribution of energy. The discrep- 
ancy in redistribution arises from different planet-star 
contrast ratios. The 8 fj,m IRAC point of Charbonneau 
et al (2008) is > 2a different from the IRS spectrum, 
and from the 8 /im IRAC measurement of Knutson et al 
(2007). Similar to the IRS spectrum, the Spitzer broad 
band photometry shows no sign of a thermal inversion. 

The Spitzer broadband photometry enables significant 
constraints on molecular abundances. Again, we take the 
£ 2 = 2 surface, corresponding to an average difference of 
1.4 a between the model and data. We find the lower 
limit of the H 2 mixing ratio to be 10~ 6 , and the upper 
limit of the CH4 mixing ratio to be 10 -5 (Figure [3}. 
At the £ 2 = 2 surface there are no constraints on the 
presence or abundances of CO and CO2. At the £ 2 = 
1 surface, the best fit models constrain the molecular 
mixing ratios as follows: 10~ 5 < H 2 < 10~ 3 ; CH 4 < 



2xl0~ 6 ; 7xl0~ 7 < C0 2 < 7xl0~ 5 ; CO is unconstrained 
because it is degenerate with CO2. The identification of 
CO2, at the £ 2 = 1 surface, is intriguing, and this is the 
first time CO2 is identified with Spitzer photometry. Our 
model also finds the C/O ratio to be 7 x 10~ 3 < C/O < 1 
for both the £ 2 = 1 and £ 2 = 2 surfaces. 

HST/NICMOS Spectro-photometry: The 
ffST/NICMOS data (Swain et al. 2009a) are use- 
ful because the wavelength range includes significant 
features of H 2 0, CH 4 , CO, and C0 2 (Figure©. We can 
therefore constrain the abundances even at the £ 2 = 2 
surface (Figure HJ). We find mixing ratios of: H 2 
- 10~ 4 ; CH 4 < 6 x 10~ 6 ; 2 x 10~ 4 < CO < 2 x 10~ 2 ; 
and CO2 ~ 7x 10~ 4 . It is remarkable that at the £ 2 = 2 
surface we can constrain the abundances so stringently. 
The exception is CH4, where the upper limit indicates 
that methane is not present in large quantities on the 
daysidc of HD 189733b. The C/O ratio is between 0.5 
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Fig. 4. — Constraints on the daysidc atmosphere of HD 189733b. The contours show minimum § 2 surfaces in the pa ramet er space; at 
each point in a two-parameter space, the color shows the minimum possible £ 2 out of a degenerate set of solutions (see S l4.1.1l for details). 
The filled color surfaces show the constraints due to the Spitzer IRS spectrum (Grillmair ct al. 2008). The red, orange, green and blue 
colors correspond to £ 2 ranges of 1.0 - 2.0, 2.0 - 3.0, 3.0 - 4.0, and 4.0 - 5.0, respectively. The solid lines show the constraints from Spitzer 
broadband photometry (Charbonneau ct al. 2008); the black (purple) solid lines correspond to the § 2 < 2 (§ 2 < 1) surface. The dashed 
lines show the £ 2 < 2 constraints from HST NICMOS Spectro-photometry (Swain et al. 2009a). The results follow from ~9x 10 6 models. 
The mixing ratios of species are shown as ratios by number. < dT/d(logP) > is the average temperature gradient in Layer 2 of the P-T 
profile, which contributes most to the emergent spectrum. 
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Fig. 5. — Pressure- Temperature profiles explored by models for 
the secondary eclipse spectrum of HD 189733b. Since each P — T 
profile can have multiple values of £ 2 , corresponding to the several 
possible molecular compositions, we color-code each P — T profile 
with only the minimum possible value of £ 2 (as also described in 
Figure [4j. The red, orange, green and blue colors correspond to 
models that fit the IRS spectrum (Grillmair et al. 2008) with 
minimum £ 2 in the ranges of 1.0 - 2.0, 2.0 - 3.0, 3.0 - 4.0, and 
4.0 - 5.0, respectively. The purple and gray profiles correspond to 
models that fit the Spitzer broadband photometry (Charbonneau 
ct al. 2008) with minimum £ 2 in the ranges 0.0 - 1.0 and 1.0 - 
2.0 respectively, and the brown profiles show models that fit the 
NICMOS data (Swain et al. 2009a) with the minimum £ 2 in the 
1.0 - 2.0 range. 

and 1. Our inferred abundances are consistent with 
those reported in Swain et al. (2009a), with the major 
exception of CO2. While we find CO2 mixing ratios of 
7 x 10~ 4 , Swain et al. (2009a) find 10" 7 < C0 2 < 10~ 6 . 
At present we have no explanation for this discrepancy. 

A mixing ratio of 7 x 10~ 4 is seemingly high from a 
theoretical standpoint; photochemistry which is thought 
to be the primary source of CO2 in hot Jupiters allows 
mixing ratios of CO2 up to ~ 10~ 6 , assuming solar metal- 
licities (Liang et al. 2003, Zahnle et al. 2009). One pos- 
sible explanation of a large inferred CO2 mixing might 
be a high mixing ratio concentrated in only a few layers 
of the atmosphere. The data can be fit equally well with 
the same mixing ratio (7 x 10~ 4 ) of CO2 present only 
between 10~ 3 - 0.3 bar causing an average mixing ratio 
of 10 -6 over the entire atmosphere. 

4.1.2. Pressure-Temperature Profiles 

Our results include constraints on the ID averaged P- 
T profile. We explore the parameter space with a grid of 
~ 10 4 P-T profiles, shown in Figure [5] For P-T profiles 
with no thermal inversions, the parameters that influence 
the spectrum most are the two pressures (Pi and P3), 
which contribute the temperature differential in layer 2 of 
the atmosphere, and the corresponding slope (a 2 ). Ad- 
ditionally, P3 governs the base flux, corresponding to the 
blackbody flux of temperature T3, on which the absorp- 
tion features are imprinted. 

The constraints on the parameters of the P-T profiles 
(§ 12. 2[) from the three data sets are shown in Figure 
The IRS spectrum (Grillmair et al. 2008) constrains the 
P-T profiles towards larger values of 02, leading to less 
steep temperature gradients. For the same base temper- 



ature (T3) of the planet dayside, the maximum tempera- 
ture gradients allowed by the broadband photometry and 
the HST spectrum are always higher than that allowed 
by the IRS data set. The £ 2 — 2 surface for the IRS data 
set constrains c*2 to be greater than 0.2, while the con- 
straints due to the broadband photometry and the HST 
data are 0.18 - 0.35 and 0.15 - 0.27, respectively. The 
same surface constrains the mean temperature gradient 
(< dT/dlogP >) in layer 2 to 100 - 400 K, 180 - 440 K 
and 290 - 470 K for the IRS, the broadband photometry 
and the HST data sets, respectively. The £ 2 = 2 sur- 
face also shows that the IRS data allows for lower base 
temperatures as compared to the other data sets. The 
skin temperature (To) at the top of the atmosphere (at 
P ~ 10~ 5 bar) is relatively unconstrained because of its 
degeneracy with the temperature gradient, and because 
of the low contribution of spectral features at that pres- 
sure. 



4.1.3. Albedo and Energy Redistribution 

Our results also constrain the albedo and day-night 
energy redistribution. The relevant parameter is 77 = 
(1 — Ab)(1 — f r ), where As is the Bond albedo, and f r is 
the fraction of input stellar flux redistributed to the night 
side (see § 13. 3[) . Since we can only constrain the prod- 
uct 77 and not the individual components, Ab and f r , one 
can estimate the constraint on A b given f r or vice- versa. 
This is similar to other forward models in the literature, 
where an f r is assumed, and the albedo is determined by 
the model. Figure |4] shows the £ 2 surface in the space of 
T e ff vs. 77. The IRS data set allows for efficient day-night 
circulation on this planet, with 77 in the range 0.38 - 0.99 
at the £ 2 = 2 level. While 77 = 1 means zero albedo and 
no redistribution, a value of rj = 0.38 translates into an 
upper-limit of f r < 0.62 (for zero albedo). Even if one 
were to assume Ab = 0.3 (a relatively high value for hot 
Jupiters; see Rowe et al. 2008), f r can be as high as 
0.46, indicating very efficient redistribution. The broad- 
band photometry and the HST data sets, on the other 
hand, require much lower efficiency in energy redistribu- 
tion, constraining the lower bound on 77, at the £ 2 = 2 
level, to values of 0.74 and 0.72, respectively. These val- 
ues translate into upper limits on f r of 0.26 and 0.28, 
respectively, suggesting low redistribution. The range of 
effective temperatures allowed by the different data sets 
follow accordingly. While the IRS data set allows for T e g 
between 1220 - 1550 K, the broadband photometry and 
HST data sets constrain T cS between 1440 - 1560 K. 

It is natural to ask: what can we say with the combined 
data sets from Spitzer IRS, Spitzer IRAC, and fiSTNIC- 
MOS? The answer depends on which error surface one is 
willing to adopt. Each of the data sets places a different 
constraint on the planet atmosphere (at the £ 2 < 2 sur- 
face). Thus, if we take the data sets at face value, and 
adopt the model £ 2 < 2 surfaces, we find that the atmo- 
sphere must be variable, both in its energy redistribution 
state and in the chemical abundances. Specifically, the 
high CO2 value could be pointing to local changes due 
to transient photochemical effects and atmospheric flows. 
On the other hand, inference of stellar variability, or re- 
fined observational uncertainties may potentially change 
this situation. 
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Fig. 6. — Transmission spectra for HD 189733b. The black filled circles with error bars show the data from Swain et al. 2008b (see S 14.21 
for details). The red, orange, and green spectra show ten models each in the £ 2 ranges of 1.0 - 2.0, 2.0 - 3.0, and 3.0 - 4.0, respectively. 
The blue spectrum shows one best-fit model, with the light-blue filled circles showing the corresponding model points binned to the grism 
bandpass, after convolving with the instrument point spread function as described in Swain ct al. (2008b). 



4.2. Transmission Spectra 

The transmission spectrum of HD189733b was pre- 
sented by Swain et al. (2008b), who reported the detec- 
tion of CH4 and H 2 in the 1.4 /im - 2.5 /im range cov- 
ered by the fiST/NICMOS spectrophotometry. Trans- 
mission data arc complementary to secondary eclipse 
data, because a transmission spectrum probes the part of 
the planet atmosphere in the vicinity of the limb, whereas 
secondary eclipse data probes the planetary day side. In 
addition, a transmission spectrum reveals line features of 
the species in absorption, imprinted on the stellar flux, 
with the self-emission of the planet contributing negligi- 
bly to the observed flux. 

We computed ~ 7 x 10 5 models to place constraints 
on allowed limb P-T profiles and quantify the allowed 
abundances by the £ 2 surfaces. In order to compute the 
£ 2 , our models were convolved with the instrument point- 
spread function, integrated over the grism transmission 
function, and finally binned to the same bins as the data 
(as described in Swain et al. 2008b). We find that 2 
of 18 binned data points (at 1.888 /im and 2.175 /im ) 
as reported by Swain et al. (2008b) have unusually high 
values, and amount to an unusually large contribution to 



the net £ 2 . In the framework of our models, we do not 
find any spectral features that might be able to explain 
those values. Consequently, we leave these two points 
out of our present analysis. Figure [6] shows the model 
spectra and data. 

4.2.1. Molecular Abundances 

Our results confirm the presence of H2O and CH4 in 
the limb, as also reported by Swain et al (2008b). At the 
£ 2 = 2 surface, as shown in Figure the H2O mixing 
ratio is constrained to 5 x 10~ 4 — 0.1, and the CH4 mix- 
ing ratio lies between 10 -5 — 0.3. These ranges of mixing 
ratios for H2O and CH4 are consistent with the mixing 
ratios found for the best-fit model reported by Swain et 
al (2008b). Although, these constraints are roughly con- 
sistent with the constraints on the same molecules on the 
dayside inferred in § 14.11 the high mixing ratios of CH4 
allowed is appaling at first glance. This hints at a region 
of high CH4 concentration in the limb, possibly on the 
night side. The effect of other possible molecules, CO, 
NH3 and CO2 on the spectrum is minimal. Although 
each of these molecules have some spectral features in 
the observed wavelength range, the upper-limits on their 
compositions placed by the data are not statistically sig- 
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Fig. 7. — Constraints on atmospheric properties at the limb of HD 189733b. The contours show mini mum £ 2 surfaces (see caption of 
Figuref4]for details), as constrained by the transmission spectrum of Swain ct al. (2008b) described in § 14.21 The red, orange, green and 
blue surfaces correspond to £ 2 in the ranges 1.0 - 2.0, 2.0 - 3.0, 3.0 - 4.0, and 4.0 - 5.0, respectively. 

the temperature T3 at the limb is constrained to lie be- 
tween 1500 K - 2250 K. The mean temperature gradient 
(< dT/dlogP >) in layer 2 constrained by the £ 2 = 2 
surface lies between 130 K - 550 K. 

5. HD 209458B RESULTS: SECONDARY ECLIPSE 

The dayside of HD 209458b has been of substantial in- 
terest, owing to the indication of a thermal inversion in 
the atmosphere of this planet. Doming ct al. (2005) 
reported the first detection of thermal emission from 
HD 209458b in the 24 (im MIPS Channel of Spitzer. 
Knutson et al. (2008) reported broadband photometry 
of the planet-star flux ratio in the four Spitzer IRAC 
channels (3.5, 4.6, 5.8, 8 //m). Newer photometric ob- 
servations have been made by Deming in the 16 fjxn IRS 
Channel and, an updated value is available for the 24 fjxa. 
MIPS Channel (Deming, private communication 2009). 
Richardson et al. (2003) observed an upper-limit on the 
contrast ratio in the K-band (2.2 ^m) using IRTF SpeX. 
Additionally, Richardson et al. (2007) and Swain et al. 
(2008a) presented a low S /N Spitzer IRS spectrum in the 
7.5-13.2 /urn range. And, after submitting our present 
work, we became aware of an independent inference of 
H 2 0, CH 4 , and C0 2 in HD 209458b, based on HSTNIC- 
MOS observations, and five-channel Spitzer photometry; 
although, a quantitative measure of fit was not reported 
(Swain et al. 2009b). Here, we report independent detec- 
tions of H2O, CH4, CO2, and CO based on quantitative 
constraints at the £ 2 < 1 level of model fits to six-channel 
Spitzer photometry. 

For HD 209458, we compute model fits to the six 
Spitzer broadband photometric detections (at 3.6, 4.5, 
5.8, 8, 16 and 24 fim), combining data from Knutson et 




Fig. 8. — Pressure- Temperature profiles explored by models for 
the transmission spectrum of HD 189733b. The profiles are color- 
coded by minimum £ 2 (see Figure [§}, and the colors are described 
in Figure [7] 



nificant. 

4.2.2. Pressure- Temperature Profiles 

Figure |8] shows the P-T profiles explored for the 
HD 189733b transmission data set. The best-fit mod- 
els with £ 2 < 2 are shown in red. Figure [7] shows the £ 2 
surface in the parameters space of the P-T profiles. The 
best-fit models arc consistent with an atmosphere with 
no thermal inversions. The ai parameter is constrained 
by the £ 2 = 2 surface to be greater than ~ 0.15, and 
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Fig. 9. — Secondary Eclipse spectra for HD 209458b. The black filled circles with error bars show the data, obtained by Spitzer photometry. 
The 3.6 fim, 4.5 fim, 5.8 /im, and 8 fim data are from Knutson ct al. (2008). The 16 /im and 24 (im data are from Doming (personal 
communication, 2009). The upper-limit at 2.2 (im shows the l-cr constraint from Richardson et al. (2003). The red, orange, and green 
spectra correspond to models with £ 2 in the range of 0.0 - 2.0, 2.0 - 3.0, and 3.0 - 4.0, respectively. For guidance, a single best fit model 
spectrum is shown in blue (with the bandpass averaged points shown in light- blue). 



al. (2008) and Deming (private communication, 2009). 
We do not use the IRS spectrum (Richardson et al. 2007; 
Swain et al. 2008a) because the S/N is too low to yield 
signficant model constraints. 

We explore a grid of 6 x 10 6 models for this system. 
Figure M shows sample model spectra and the data for 
the planet-star flux contrast in this system. The figure 
shows model spectra at different levels of £ 2 . As shown in 
the figure, our best fit spectra (a sample binned spectrum 
is shown in light blue) fit all the available data to within 
the l-cr error bars. And, our £ 2 = 2 models match the 
observations to ~ 1.4er on average. 

5.1. Molecular Abundances 

Our results indicate the presence of H 2 0, CO, CH 4 
and C0 2 in the atmosphere of HD 209458b. As will be 
discussed in § 15.21 the allowed P-T profiles show the ex- 
istence of a deep thermal inversion in the atmosphere 
of HD 209458b. Consequently, all the line features of 
the molecules are seen as emission features, rather than 
absorption features as in case of HD 189733b. The con- 
straints on the concentrations of all the molecules arc 
shown in Figure 1101 



The £ 2 = 2 surface places an upper-limit on the mixing 
ratio of H 2 at 10~ 4 , with the £ 2 = 1 surface allowing 
mixing ratios up to 10~ 5 . The limits on H 2 are gov- 
erned by the band features of the molecule in the 5.8 /jm 
and 24 fim channels. While the upper limit on the mix- 
ing ratio of H 2 is well constrained, the constraint on 
the lower limit is rather weak, allowing values as low as 
10 -8 at £ 2 = 1. The stringent upper limit on H 2 comes 
primarily from the low planct-star flux contrast observed 
in the 24 /im channel. Similarly, the high flux contrast 
observed in the 5.8 fim channel should, in principle, yield 
a stringent lower limit. However, the weak lower limit is 
because of the large observational uncertainty on the flux 
contrast in the 5.8 /im channel. The 1 — a error bar on 
the observed 5.8 /jm flux contrast is the largest of uncer- 
tainties in all the channels, and 1.5 times larger than the 
error bar on the 24 fim flux contrast. 

The £ 2 = 2 constraint on the mixing ratio of CH4 is 
10~ 8 - 0.04, based on band features in the 3.6 /im and 8 
fim channels. Although, the features in the 3.6 fim chan- 
nel are degenerate with some features of H 2 in the same 
channel, the exclusive CH4 features in the 8 fim channel 
break the degeneracy. The £ 2 = 2 surface places an abso- 
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Fig. 10. — Constraints on the dayside atmosphere of HD 209458b. The contours show minimum § 2 surfaces (see Figure[4j in the parameter 
space of models for secondary eclipse spectra of HD 209458b. The red, orange, green and blue surfaces correspond to § 2 in the ranges 1.0 
- 2.0, 2.0 - 3.0, 3.0 - 4.0, and 4.0 - 5.0, respectively. Additionally, the purple colored surfaces show regions that allow £ 2 < 1. The results 
follow from ~6x 10 6 models with thermal inversions. Models without thermal inversions are ruled out by the data. < dT/d(logP) > is the 
average temperature gradient in the atmosphere between pressures Pi and P2, in P-T profiles with thermal inversions, which contributes 
most to the observed emission features. 
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Fig. 11. — Pressure- Temperature profiles explored by models for 
the secondary eclipse spectrum of HD 209458b. All the profiles 
have an inversion layer. Profiles without thermal inversions are 
ruled out by the data. The profiles are color-coded by minimu m 
(as described in Figurc[5]l, and the colors are described in Figurc[TO] 

lute lower-limit on CO at 4x 10~ 5 , and requires a mixing 
ratio of C0 2 between 2 x 1CT 9 - 7 x 1CT 6 . The need for 
CO2 arises from the high contrast observed between the 
4.5 //m and the 3.6 /im channels. The 4.5 /jm channel 
could in principle be explained by CO alone. However, 
such a proposition would require large amounts of CO, 
close to mixing ratios of 1. CO2 provides additional ab- 
sorption in the same 4.5 /jm channel, with a relatively 
reasonable abundance. In addition, the amount of CO2 
is also constrained by the 16 /jm channel — too high 
of a CO2 abundance would not agree with the data be- 
cause C0 2 has a strong absorption feature in the 16 fim 
channel. Thus, the 4.5 /j,m and 16 fim channels together 
constrain the amounts of CO2 , thereby also constraining 
the CO feature in the 4.5 (im channel. If the 16 /im data 
did not exist, there would have been a large degeneracy 
in the abundances of CO and CO2. 

If we consider the £ 2 = 1 surface, instead of the £ 2 = 2 
surface, our results place tighter constraints on all the 
molecules. At this surface, the constraints on the mixing 
ratios are: 1CT 8 < H 2 < 1CT 5 ; 4 x 1CT 8 < CH 4 < 0.03; 
CO > 4 x 10~ 4 ; and, 4 x 10~ 9 < C0 2 < 7 x 10~ 8 . 
However, in considering the £ 2 = 1 surface, we might be 
over-interpreting the data. 

We find definitive detections of CH4, CO and C0 2 . 
Even at the £ 2 = 4 surface, the data requires non-zero 
mixing ratios for the carbon-bearing molecules. Further- 
more, the results suggest a high C/O ratio (see Seager et 
al. 2005). Although the model fits give C/O ratio close 
to 1 or higher, our molecular abundance grid resolution 
is not high enough to put a precise lower limit to the C/O 
ratio. As seen from Figure [THl the allowed models reach 
values as high as ~ 60 for £ 2 = 1, and ~ 600 for £ 2 = 2. 
As discussed previously and in § 6.1, high molecular con- 
centrations and high C/O ratios could be localized effects 
overrepresented by the ID averaged retrieval. 

5.2. Temperature Structure 

The most interesting aspect about the atmosphere of 
HD 209458b is the presence of a thermal inversion on the 



planet dayside, as previous studies have pointed out (e.g., 
Burrows et al. 2008). While modeling atmospheric spec- 
tra, one often encounters a degeneracy between temper- 
ature structure and molecular compositions. However, 
by exploring a large number of P-T profiles we find that 
the day-side observations for HD 209458b cannot be ex- 
plained without a thermal inversion, for any chemical 
composition, at the xi 2 < 2 level. On the other hand, 
our results place stringent constraints on the ID averaged 
structure of the required thermal inversion. 

Our results presented here explore ~ 7000 P-T pro- 
files with thermal inversions. Figure [TT] shows all the 
P-T profiles with £ 2 < 5. The ranges in the P-T param- 
eters explored are shown in Figure [TUJ In the presence 
of a thermal inversion, the spectra are dominated by the 
location and extent of the inversion layer, quantified by 
the parameters P2, T 2 , and a.2- We find the models with 
£ 2 < 2 to have an inversion layer in the lower atmosphere, 
at pressures (P2) above ~ 0.03 bar, with the best-fit mod- 
els preferring the deeper layers, i.e., higher P 2 - We find 
that the £ 2 < 1 region extends to inversion layers as deep 
as P2 ~ 1 bar. However, we limit P2 to a value of 0.5 
based on the physical consideration that, at pressure > 
1 bar collisional opacities are likely to dominate over the 
molecular opacities in the optical which are expected to 
be the cause of thermal inversions. And, therefore, we 
do not expect inversions at such high pressures. 

The allowed models (Figure IT0|) cover a wide range 
in the base temperature (T3) owing to the fact that the 
spectrum is dominated by the maximum temperature dif- 
ferential (Ti - T2) attained in the stratosphere, and is less 
sensitive to T3 by itself. The stratospheric temperatures 
of all the models with £ 2 < 2 lie between 800 K and 
2200 K. The temperature gradient in the inversion layer 
is governed by the U2 parameter. As can be seen from 
Figure[Tni the data restricts the £ 2 < 2 models to have 0:2 
between 0.05 and 0.25, with the £ 2 < 1 models preferring 
values between 0.05 and 0.15. These values of ct2 trans- 
late into mean temperature gradients (< dT / d log P >) 
in the inversion layer to be 230 K - 1480 K for £ 2 < 2, and 
400 K - 1480 K for £ 2 < 1. Consequently, large thermal 
inversions spanning more than a few dex in pressure, and 
hence with less steeper thermal gradients, arc ruled out 
by the data. This is also evident from Figure QTJ Finally, 
£ 2 is not very sensitive to the parameters in Layer 1 of 
the P-T profile, since the dominant spectral features in 
this planet arise from the inversion layer. 

5.3. Albedo and Energy Redistribution 

Our results also show the constraints on the albedo and 
day-night redistribution. Figure [TU] shows the £ 2 surface 
in the space of T e ff vs. 77, where, 77 = (1 — Ab){1 — f r )- 
As is the Bond albedo, and f r is the fraction of input 
stellar flux redistributed to the night side. As described 
in § I4.1.3[ using 77 one can estimate the constraint on 
A B given f r , or vice- versa. For HD 209458b, the £ 2 = 2 
surface in Figure flUl constrains rj in the range 0.24 - 0.81, 
with the £ 2 = 1 surface in the range 0.26 - 0.72. If we 
assume an A b of 0.24 (2-cr upper-limit, scaled from the 
1-er estimate of Rowe et al. 2008), the constraints on f r 
for the £ 2 = 2 and £ 2 = 1 surfaces are 0.0 - 0.68 and 0.05 
- 0.66, respectively. This range of values allows for very 
efficient redistribution in the atmosphere of HD 209458b, 
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although values of f r > 0.5 seem physically unrealistic. 

6. DISCUSSION 

Our results in this work have been presented as 
goodness-of-fit contours in the parameter space of a ID 
averaged atmosphere model. In this section, we first dis- 
cuss the interpretation of our results in the light of the 
fact that the atmosphere is inherently 3D, with dynam- 
ical effects that cannot be captured in a ID model. We 
then present an assessment of our best-fit models vis- 
a-vis the goodness-of-fit contours. Finally, we discuss 
possible future extensions to our model. 

6.1. Contribution Functions 

We have constrained the P-T profile and abundances 
on the planetary day side. But which part of the at- 
mosphere is really being sampled by the observations? 
Here we will discuss the contribution to the net emer- 
gent flux from the surface elements at the top of the 
atmosphere (see Knutson et al. 2009b for a discussion of 
radial contribution functions). We begin with the sim- 
ple picture of radiative equilibrium, with no atmospheric 
circulation, in a homogeneous and isotropically radiating 
atmosphere. In this case, the amount of flux from a given 
surface clement is proportional to the amount of incident 
flux from the star, where the incident stellar flux drops 
off as cos 0; 8 is the angle away from the sub-stellar point. 
In this picture, the planetary infrared flux is greatest at 
the sub-stellar point, where the planet is being heated 
the most. The area contribution to the net planet flux, 
however, increases as smO. The bigger annuli correspond 
to larger areas on the projected planetary disk as viewed 
by the observer. The lower flux emitting regions (the 
outer regions of the disk, with lower incident stellar in- 
tensity), therefore, have more area than the higher flux 
regions. There is a sweet spot at 45 degrees where the 
contribution to the total observed flux peaks, with a dis- 
tribution of sin(26>). In other words, the total flux is an 
average over a wide region on the projected planetary 
disk. 

A more realistic scenario is one not limited to radia- 
tive equilibrium; indeed atmospheric circulation models 
show dramatic differences of emergent flux across the 
planetary disk caused by strong hydrodynamic flows (see, 
e.g., Showman et al. 2009; Langton and Laughlin, 2008; 
Cho et al. 2003; Dobbs-Dixon and Lin. 2008). Without 
knowing the atmospheric circulation and the resulting 
surface pattern, can we say anything about which part of 
the day-side planetary disk our P-T profiles and abun- 
dances preferentially sample? Given a retrieved abun- 
dance for a molecule, we do not know if this is from a 
localized concentration or from a global abundance, or 
something in between. For instance, we have found ev- 
idence for a high CO2 abundance in HD 189733b (see 
also Swain et al. 2009a). But is this just a local effect 
of transported photochemical products? Or is the whole 
day-side evenly enriched in CO2? If we are not measur- 
ing a true average, we have to be careful about how we 
interpret the retrieved CO2 abundance. 

A thorough examination of the contribution of differ- 
ent parts of the planetary disk has to take into account 
the emergent surface flux profiles that are more realis- 
tic than the assumption of homogeneity and radiative 



equilibrium. For example, 3-D models of hot Jupiter at- 
mospheres by Showman et al. (2009) show a much slower 
than cos 6* drop off from the sub-stellar point caused by 
hydrodynamic flows. And at the planetary photosphere 
the "sub-stellar hot spot" is blown downwind, to 30 de- 
grees; a detail that would not be apparent in the as- 
sumption of an isotropic and homogenous atmosphere. 
Moreover, the Spitzer 8 micron measurements by Knut- 
son et al. (2007) show a hotspot at lower separation from 
the sub-stellar point (and also show a cold spot on the 
same hemisphere). For instance, if CO2 were predomi- 
nantly formed in localized regions and blown downwind, 
it could land at geometrically favored regions in terms of 
the ID averaged spectrum. Alternatively, vertically lo- 
calized concentrations of CO2 at favorable pressure levels 
could give rise to strong spectral signatures, similar to 
those formed by having the same concentration of CO2 
over the entire atmosphere. In these cases, a high value 
of CO2 may not be representative of the entire day side. 

6.2. Best Fit Models vs. Error Surface 

In the current work we have described most of our 
results at the £ 2 = 2 level. Our formulation for £ 2 , as de- 
scribed in § 13.51 would be equivalent to the conventional 
definition of reduced \ 2 1 if ^obs were to be replaced by 
the number of degrees of freedom. Normally, one aims 
for fits at the reduced x 2 = 1 level (we have also provided 
results for fits at the level of £ 2 = 1). A \ 2 , however, is 
devoid of meaning when the number of degrees of free- 
dom is less than zero — in that we cannot relate x 2 to a 
confidence level. In the absence of confidence levels, we 
are considering £ 2 as a weighted mean squared error. In 
this framework, a £ 2 = 1 means that, on average, the 
model deviates from the mean values of the data by 1-a, 
and £ 2 = 2 means that the deviation is \Jl = 1.4 o. In 
this work, the number of data points in the Spitzer broad- 
band photometry is less than the number of free model 
parameters, meaning less than zero degrees of freedom. 
In both the Spitzer IRS and the fiST/NICMOS data sets, 
there are, in principle, more data points than the num- 
ber of free model parameters. But both these data sets 
sample a narrow wavelength range, with a limited num- 
ber of absorption features. In addition, the degeneracies 
in our model parameter space are poorly understood, if 
at all. Without a robust parameter-fitting algorithm, we 
cannot assign confidence levels even with the HST and 
IRS data sets. Which error surface should be used to 
encompass the best fit models, given that we have no 
direct correspondence to confidence levels? There is no 
conclusive answer and, ideally, one might like to consider 
only those models which fit the data to within £ 2 = 1, 
or even £ 2 = indicating perfect fits. However, we have 
used a conservative choice of the £ 2 = 2 surface for inter- 
pretation of our results, allowing for models fitting the 
observations to within ~ 1.4(7 of each data point, on an 
average. 

One approach to accurately consider the uncertainty 
on the data is as follows. One could take tens of thou- 
sands of realizations of each data set, uniformly sampled 
over the error distribution. Fitting each realized data 
set would mean running millions of models for each case, 
leading up to billions of models per data set per planet. 
Such an approach is computationally prohibitive at the 
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present time. 

6.3. Suggested Variability 

There is evidence that the dayside atmosphere of 
HD 189733b is variable. One highly relevant example 
to our interpretation is the 8 [im Spitzer IRAC measure- 
ment of HD 189733b reported by Knutson et al. (2007) 
versus that reported by Charbonneau et al. (2008) in the 
same channel at a different time. The flux ratios differ 
by ~ 2.3 <7. Also significant is the disagreement in the 
planet-star contrast levels between the HD 189733b IRS 
spectrum and the IRAC broadband photometry at 8 /im, 
but not at 5.8 ^m (Grillmair et al. 2008). 

Atmospheric variability is also suggested by the very 
strong C0 2 feature in the HD 189733 HST/NICMOS 
data set and the lack of a noticeably strong absorption 
feature in the 4.5 /xm and 16 /jm Spitzer channels. It 
is well known that the CO2 molecule has an extremely 
strong absorption feature at 15 /im; yet there is no indi- 
cation of strong CO2 absorption observed in the 16 fj,m 
channel. Whether this data variability is due to intrin- 
sic variability of the planet atmosphere or instead is due 
to the different data reduction techniques is under study 
(Desert et al. 2009; Beaulieu et al. 2009). 

Turning to our model results, we have been unable to 
find fits to all three data sets of HD 189733b at the £ 2 < 2 
level, in the large range of parameter space we have ex- 
plored. There may be befitting P-T profiles and com- 
positions that our grid did not cover, but we think this 
is unlikely, mainly because the data sets show noticeable 
differences even without running the models, e.g. the 
observed features of CO2. We plan to implement a more 
exhaustive exploration of the P-T parameter space by 
extending our work with a new optimization technique 
in the future. 

6.4. Model Extensions 

Our approach allowed us to run a large ensemble of 
models in the parameter space. Nevertheless, future cal- 
culations are needed to estimate the covariances between 
the model parameters. Several model parameters are de- 
generate; for example, molecular abundances of species 
that share absorption features in the same photometric 
channel (e.g., CO and CO2). Also well-known is the 
degeneracy between the P-T profile and the molecular 
abundances. Understanding the covariances between the 
parameters might help break the degeneracies in the fu- 
ture. 

We chose the ~ 10 4 P-T profiles based on physical ar- 
guments and preliminary fits to the data. Even though 
having a grid of ~ 10 4 P-T profiles in the six-parameter 
P-T space, and ~ 10 7 models overall, means that our 
grid is coarse, we believe it covers most physically plau- 
sible profiles. Nonetheless, a more complete exploration 
requires a more sophisticated approach, and may find P- 
T profiles and the corresponding abundance constraints 
that may be slightly different from those reported here. 

Our model does not include clouds or scattering. Al- 
though cloud-less models are justified for hot Jupitcrs 
(see § 13. 2[) . clouds are needed for extending our model 
to cooler planetary atmospheres, such as those of hab- 
itable zone super Earths and cooler giant planets. Our 
model also does not include photochemistry, but instead 



varies over different molecular abundances. Once a set 
of best-fit models are found by our approach, the molec- 
ular abundances can be used to guide photochemical 
models; conversely, photochemical models can be used 
to check the plausibility of the best-fit molecular abun- 
dances. And, although we have used the latest avail- 
able high-temperature opacities, further improvements 
are sought, particularly for CH4 and CO2. 

Our model presented in this work is suitable only for in- 
terpreting the infrared spectra of planetary atmospheres; 
it was motivated by existing observations. We do not 
include the sources of opacity and scattering that are re- 
quired to model the visible part of the spectra. However, 
our model can be extended to visible wavelengths, and 
we aim to discuss it in a future work. Constraints placed 
by visible observations will help to further narrow down 
the model parameter space already constrained by IR 
observations. 

We have not included the intrinsic flux from the planet 
interior, because, for hot Jupiters stellar irradiation is 
the dominant energy source. However, extensions of our 
model to brown dwarf atmospheres or young, or mas- 
sive, Jupiters will require inclusion of the interior en- 
ergy. And, although irradiation is naturally included in 
the form of the parametric P-T profile, our model gives 
no information about the nature of the absorber. Never- 
theless, given the best-fit P-T profile from our approach, 
one could quantify the range of parameters of a potential 
absorber. 

6.5. Temperature Retrieval as a Starting Point for 
Forward Models 

We have presented a radical departure from the stan- 
dard exoplanet atmosphere models. Our temperature 
and abundance retrieval method does not use the as- 
sumption of radiative + convective equilibrium in each 
layer to determine the P-T profile. By searching over 
millions of models, we are extracting the P-T profile from 
the data. 

We have not necessarily suggested to replace so-called 
"self-consistent forward models", those models that in- 
clude: radiative transfer, hydrostatic equilibrium, radia- 
tive + convective equilibrium, day-night energy redistri- 
bution, non-equilibrium chemistry, and cloud formation. 
Instead of a stand-alone model, our method can be used 
as a starting point for narrowing the potential parameter 
space to run forward models in. We anticipate that this 
temperature retrieval method will enable modelers to run 
enough forward models to find quantitative constraints 
from the data, instead of running only a few representa- 
tive models. 

Atmospheric temperature and abundance retrieval is 
not new. Solar system planet atmosphere studies use 
temperature retrieval methods, but in the present con- 
text there is one major difference. Exquisite data for so- 
lar system planet atmospheres means that a fiducial P-T 
profile can be derived. The temperature retrieval process 
for a solar system planet atmosphere therefore involves 
perturbing the fiducial profile. For exoplanets, there is 
no starting point to derive a fiducial model because of the 
limited data. Building on methods used for solar system 
planets, Swain et al. (2008b & 2009) have used a temper- 
ature retrieval method using published P-T profiles from 
other models, and perturbing the profiles to find fits to 
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data (also see Sing et al. 2008). In this work, we report 
a new approach for atmospheric retrieval of cxoplancts 
which allows us to explore a wide region of the parame- 
ter space, and to possibly motivate a new generation of 
retrieval methods tailored to exoplanet atmospheres. 

7. SUMMARY AND CONCLUSIONS 

We have presented a solution to the retrieval prob- 
lem in exoplanet atmospheres: given a set of obser- 
vations, what are the constraints on the atmospheric 
properties? Our method consists of a new model in- 
volving line-by-line radiative transfer with a paramet- 
ric pressure-temperaure (P-T) profile, along with physi- 
cal constraints of hydrostatic equilibrium and global en- 
ergy balance, and allowing for non-equilibrium molecular 
compositions. This approach enabled us to run millions 
of models of atmospheric spectra over a grid in the pa- 
rameter space, and to calculate goodncss-of-fit contours 
for a given data set. Thus, given a set of observations and 
a desired measure of fit, our method places constraints 
on the P-T profile and the molecular abundances in the 
atmosphere. 

We used our method to place constraints on the at- 
mospheric properties of two hot Jupiters with the best 
available data: HD 189733b and HD 209458b. We com- 
puted ~ 10 7 models for each system. Our constraints on 
the P-T profiles confirm the presence of a thermal inver- 
sion on the dayside of HD 209458b, and the absence of 
one in HD 189733b, as has been previously reported in 
the literature. And, the constraints on molecular abun- 
dances confirm the presence of H 2 0, CH4, CO and CO2 
in HD 189733b, as has been found by several other stud- 
ies in the past. 

We presented some new findings. We report indepen- 
dent detections of H2O, CH4, CO2, and CO, from six- 
channel Spitzer photometry of HD 209458b, at the £ 2 < 1 
level (also see Swain et al. 2009b). At the £ 2 = 2 surface, 
we find the lower-limits on the mixing ratios of CH4, CO 
and C0 2 to be 10" 8 , 4 x 10~ 5 and 2 x 10~ 9 , respectively. 
In addition, we find an upper-limit on H2O mixing ra- 
tio of 10~ 4 . Secondly, our results using the high S/N 
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IRS spectrum of HD 189733b indicate that the best fit 
models are consistent with efficient day-night circulation 
on this planet, consistent with the findings of Knutson 
et al. (2007 & 2009b); although, we also find that the 
broad band photometry data requires much less efficient 
day-night circulation, as has been reported by several 
studies in the past. The different constraints placed by 
the different data sets arise from the variability in the 
data. 

A few key observations of HD 189733b vary between 
different data sets at similar wavelengths. Moreover, a 
noticeably strong CO2 absorption in one data set is sig- 
nificantly weaker in another. We must, therefore, ac- 
knowledge the strong possibility that the atmosphere is 
variable, both in its energy redistribution state and in 
the chemical abundances. 

Our new temperature and abundance retrieval — to- 
gether with the limited but remarkable data sets - 
charts a new course for exoplanet atmosphere interpre- 
tation. Whether our approach is used as a stand alone 
model, or as a starting point for forward models, the field 
of exoplanet atmospheres is ready for a technique to place 
detailed atmospheric constraints governed by the data. 
The ideal data set is one with a wide wavelength range 
that covers more than one spectral feature of different 
molecules. The James Webb Space Telescope, scheduled 
for launch in 2013, is ideal for observations of transiting 
exoplanet atmospheres and should usher a new era for 
quantitative analyses of exoplanet atmospheres. 
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